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Abstract 

In this paper we describe a general, computationally feasible strategy to 
deduce a family of interpolatory non-stationary subdivision schemes from 
a symmetric non-stationary, non-interpolatory one satisfying quite mild as- 
sumptions. To achieve this result we extend our previous work [C. Conti, 
L. Gemignani, L. Romani, Linear Algebra Appl. 431 (2009), no. 10, 1971- 
1987] to full generality by removing additional assumptions on the input sym- 
bols. For the so obtained interpolatory schemes we prove that they are capable 
of reproducing the same exponential polynomial space as the one generated 
by the original approximating scheme. Moreover, we specialize the compu- 
tational methods for the case of symbols obtained by shifted non-stationary 
affine combinations of exponential B-splines, that are at the basis of most non- 
stationary subdivision schemes. In this case we find that the associated family 
of interpolatory symbols can be determined to satisfy a suitable set of general- 
ized interpolating conditions at the set of the zeros (with reversed signs) of the 
input symbol. Finally, we discuss some computational examples by showing 
that the proposed approach can yield novel smooth non-stationary interpola- 
tory subdivision schemes possessing very interesting reproduction properties. 
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1 Introduction 



Binary interpolatory subdivision schemes are efficient iterative procedures for the 
generation of interpolatory curves: starting with the set of points to be interpolated, 
at each recursion step a new point is inserted in between any two given points so 
that the limit curve, whenever exists, not only interpolates the initial set of points 
but also all the points generated through the whole process. Taking into account 
that a curve is displayed on the screen by visualizing a discrete set of its points, 
from a computational viewpoint interpolatory subdivision schemes turn out to be 
more efficient than classical interpolating methods in several situations. In fact the 
limit points obtained within five or six subdivision iterations are in general enough 
for a good discrete representation of the limit shape. This is one of the reasons 
why interpolatory subdivision schemes are widely used in applications and often 
preferred to standard methods. 

Two important areas where interpolatory subdivision schemes play a crucial role 
are Computer Aided Geometric Design (CAGD) and wavelets construction (see [10] 
and [18], respectively). In these fields a fundamental issue that recently emerged 
is concerned with the study of numerical algorithms for converting known approxi- 
mating schemes into new interpolatory ones. Starting from the works [T7] and [22] . 
where the conversion is obtained for a specific approximating scheme by means of a 
push-back or a tweak operator, geometric approaches based on the idea that an in- 
terpolatory refinement can be interpreted as an averaging step on the control points 
followed by a further adjustment of some of them to ffi the interpolation constraints 
were presented [151 CB] • Very recently a completely different technique relying upon 
the interplay between polynomial and structured matrix computations has been pro- 
posed in [5] . In that work for a given symmetric Hurwitz approximating symbol an 
associated family of interpolatory symbols is determined in such a way to satisfy an 
auxiliary polynomial equation. As it clearly appears, although the latter strategy 
turns out to be more general than the previous ones, it is limited to the context 
of stationary subdivision schemes. Being non-stationary subdivision schemes more 
powerful than stationary ones and very attractive in several applications such as in 
CAGD (because of their ability to reproduce conic sections, spirals or widely used 
trigonometric curves) it is of fundamental importance to provide a general and ef- 
ficient method to convert a given non- stationary, non-interpolatory scheme into a 
family of interpolatory ones. To our knowledge, there exists only a new paper [T] 
addressing this problem, which presents a strategy that is restricted to the case of 
symmetric subdivision masks of odd width, namely symmetric subdivision symbols 
of even degree. 

The goal of this paper is to elaborate on our recent work [5j to progress along 
different directions. In particular, (i) we extend the applicability of the proposed 
construction, (ii) we investigate the reproduction properties of the so-obtained in- 
terpolatory schemes and (iii) we design algorithms specifically suited for the case of 
approximating symbols generated from exponential B-splines, that are at the basis of 
most non-stationary subdivision schemes. More specifically, in this paper we prove 
that the strategy described in [S] can still be pursued under very relaxed conditions 
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on the approximating symbols we deal with, say {a^''\z), k > 0}. If, for a given 
fixed A; > 0, a^'^^z), a^^\—z) are relatively prime polynomials, then a double family 
of interpolatory symbols associated with a^'^\z) can be generated by solving two 
different Bezout-like polynomial equations. In the symmetric case where a^^^ [z) is a 
symmetric polynomial, the double family reduces to one single family since the solu- 
tions of these two equations are suitably related. In the Hurwitz case where a^^"^ [z] 
is a Hurwitz polynomial, the distribution of the roots implies the primality condi- 
tion. Whenever such a condition is satisfied for any /c > then the correspondence 
of a'^'^\z) with any member of the associated double family allows one to define a 
family of interpolatory subdivision schemes derived from the given non- stationary 
approximating one. The computation of the interpolatory symbol amounts to solve 
the corresponding polynomial equation. If the approximating symbol is specified by 
spectral information, as it is generally the case of exponential B-splines, then it is 
shown that the equation can be efficiently solved by using the tool of (incomplete) 
partial fraction decomposition. This gives a representation of the associated interpo- 
latory symbol in terms of a set of generalized interpolating conditions attained at the 
zeros (with reversed signs) of the approximating symbol. For the newly generated 
interpolating schemes we prove an important reproduction result: the exponential 
polynomial space reproduced by the interpolatory scheme is the same function space 
generated by the approximating one it is originated from. On the contrary, a general 
result concerning convergence and/or smoothness of a non- stationary interpolatory 
subdivision scheme induced by a non-stationary approximating one is not yet avail- 
able. However, in many specific examples we have considered, the analysis can be 
performed by using ad-hoc techniques. In this way, by starting with approximating 
schemes suitably generated by five term affine combinations of exponential B-splines, 
we are able to find novel smooth non-stationary interpolatory subdivision schemes 
possessing very interesting reproduction properties. 

The paper is organized as follows. In Section [2] the needed background on non- 
stationary subdivision schemes is given. In Subsection 13.11 we review and generalize 
the basic strategy proposed in [3] for the construction of an interpolatory subdi- 
vision mask from a given approximating one. Effective computational procedures 
for implementing this strategy are discussed in Subsection 13.21 These procedures 
are the key ingredients of our algorithm, named Appint and stated in Subsection 
13. 3[ to move from a non- stationary approximating subdivision scheme to a family 
of non-stationary interpolatory ones. The reproduction properties of these schemes 
are studied in Section H] whereas in Section [5] the application of the algorithm to 
several instances of non-stationary approximating subdivision schemes generating 
exponential polynomials is considered. Finally, conclusions and further work are 
drawn in Section [61 
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2 Background 



In this section we briefly recall some needed background on stationary and non- 
stationary subdivision schemes. For more material on subdivision schemes we refer 
the reader to the seminal work by Cavaretta, Dahmen and Micchelli [1] , to the more 
recent survey by Dyn and Levin [10] and to the well-known book by Warren and 
Weimer pi] . 

Subdivision schemes are simple iterative algorithms to efficiently generate curves 
and surfaces. Any subdivision scheme is defined by an infinite sequence of coefficients 
collected in the so called refinement masks k > 0}. We assume that any mask 

a^'^) := (^af'^ G M, z G is of real numbers and has finite support for all /c > i.e. 

af'^ = for i ^ [—n{k),n{k)] for suitable n{k) > 0. The k-level subdivision operator 
associated with the /c-level mask a^'^^ is 

: £(Z) ^ £(Z) , q), := Y^af'J,^ g„ ^ G Z , (2.1) 

where £(Z) denotes the linear space of real sequences indexed by Z whose elements 
will be denoted by boldface letter, q := (gj G M, i G Z). The subdivision scheme con- 
sists of the subsequent application of S'a(o), ■ ■ ■ , S'^c*) from a given starting sequence, 
say q, generating the scalar sequences 

q(0) := q , qC^^+i) ;= S^i,, q^'^^ for k > 0. (2.2) 

In case the masks {a^^^ k >0} are kept fixed over the iterations, that is a'^'^^ = a for 
all k > 0, the subdivision scheme is said to be stationary, otherwise non-stationary. 
Attaching the data generated at the k-th step to the parameter values t^''^ with 

tf)<4a, and e-f^ = 2-^ k>0 

(these are usually set as tf'^ := ^) we see that the subdivision process generates 
denser and denser sequences of data so that a notion of convergence can be estab- 
lished by taking into account the piecewise linear function Q*-'^-* that interpolates the 
data, namely 

gW(^W)^^w^ g^'^l[,(.),(.)]eni, ^GZ, k>o, 

where Hi is the space of linear polynomials. If the sequence {Q^'^\ ^ > 0} converges, 
then we denote its limit by 

U ■■= lim g(^) 

and say that /q is the limit function of the subdivision scheme based on the rule 
( 12. 2 p for the data q. Several subdivision properties can be read off from the symbols 

aW(^) = ^af^2\ A;>0, z G C \ {0} 
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associated to the masks {a^^\ k > 0}. Also, the corresponding sub-symbols 



0- 

even 



related to the symbols by the relation 

a(') {z') + z-a^''\Az^) = a^''\z), 
even^ ^ odd^ ^ ^ ^' 

are useful tools for subdivision analysis. Note that since the masks are always sup- 
posed to be finitely supported, all symbols are Laurent polynomials. Nevertheless, 
for the analysis of subdivision properties of our concern we can always assume to 
work with polynomial symbols, at least after the application of a suitable shift at 
each iteration. 

A celebrated class of stationary subdivision schemes is given by degree-n poly- 
nomial B-spline subdivision schemes, whose (unique) symbol is 

Bniz) = ^ \J , k>(}. (2.3) 

The non-stationary counterpart of (12. 3 p is the symbol of the so-called exponential 
B-splines. They are piecewise functions whose pieces are exponential polynomials 
(the latter ones will be recalled in the next definition). These are defined in terms 
of a linear differential operator and turn out to be of great interest in geometric 
modeling for the design of important analytical shapes like conic sections, spirals 
and classical trigonometric curves. 

Definition 1. (Space of exponential polynomials) LetT G Z+ and^ = (70,71, ■ ■ ■ ,7r) 
with 72- 7^ a finite set of real or imaginary numbers and let -D" the n-th order dif- 
ferentiation operator. The space of exponential polynomials Vt,7 is the subspace 

T 

Vt,j := {/ : M ^ C, / G C^{R) : J] ^jD^ f = 0}. (2.4) 

j=o 

A characterization of the space ^^,7 is provided by the following: 

Lemma 1. Let 7(2;) = X]J=o^i'^'' '^^^ denote by {^£, T^j^^i ... the set of zeros 
with multiplicity of'y{z) satisfying 

j^'\9e) = 0, r = 0,--- ,r,-l, i = l,---,N. 

It results 

N 

T = J2ri, VV,7 := 5par^{a;V^^ r = 0,--- ,r^- 1, £=1,---,N}. 
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As proved in |T9] (see also [21]) exponential B-splines can be generated via a 
non- stationary subdivision scheme based on the symbols 

B^:\z) = 2\{['-^^^;^] , k>^. (2.5) 
£=1 V + 1 / 

Its limit function belongs to the subclass of C'^"'^ degree-n L-splines [23] (with 

n = T — 1) whose pieces are exponentials of the space Vr,'y. Notice that, when = 

with Ti = n + 1, then B^n^ (z) in (12. 5p does not depend on k being the symbol of a 
degree-n B-spline given in (12. 3p . An important aspect of subdivision schemes is their 
convergence capability to specific classes of functions. In particular, a subdivision 
scheme is said to possess the property of generating exponential polynomials if, for 
any initial data uniformly sampled from some exponential polynomial function, the 
scheme yields a function belonging to the same space in the limit. Even more, the 
subdivision scheme is reproducing exponential polynomials if, for any initial data 
uniformly sampled from some exponential polynomial function, the scheme yields 
the same function in the limit. To this purpose, we recall the following two important 
definitions (see, for example, [7] and [25]). 

Definition 2 (Vy^'y- Generation). Let {a^''\z), k > 0} be a set of subdivision sym- 
bols. The subdivision scheme associated with the set of symbols {a^'^^z), k > 0} is 
said to be Vr,7-generating if it is convergent and for f e Vt,7 and for the initial 
sequence f° := {/(t°), i G Z}, it results 

lim S'a(fc) ■ ■ ■ ^^(ojf" = / , / e Vt,7 • 

k—^oo 

Definition 3 (Vr.'y-Reproduction). Let {a^^^'^z), k > 0} be a set of subdivision 
symbols. The subdivision scheme associated with the symbols {a'^{z), k > 0} is 
said to be Vr,'y-reproducing if it is convergent and for f G Vt^^ and for the initial 
sequence fO := {/(t°), i e Z}, it results 

lim Sg^(k) ■ ■ ■ Sg^(o) = f . 

Since the space of exponential polynomials trivially includes standard polynomi- 
als. Definitions [2] and [3] include, as special cases, the notion of polynomial generation 
and polynomial reproduction, respectively. For a complete analysis of the latter con- 
cepts in the stationary situation -which are very much related to the approximation 
order of the subdivision scheme- the interested reader can see jllj . 

We conclude by recalling that a subdivision scheme is said to be interpolatory if 
the refinement masks {sS''\ k > 0} satisfy 

^2? = ko, or equivalently, ^^''^^^^(2) = 1, k > 0, (2.6) 

meaning that all points generated by the subdivision process at a given level k will 
be kept in the next level k + 1. We also mention that from (12.61) it follows that a 
mask a'-'^-* is interpolatory if and only if all its symbols a^''\z) satisfy the algebraic 
condition 

a^^^ (z) + a^^^ i-z) = 2, VA; > 0. (2.7) 
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3 From approximating to interpolatory subdivi- 
sion schemes 



In this section we introduce the key ingredients of our proposed algorithm termed 
Appint to generate a family of non- stationary interpolatory subdivision schemes 
starting from an initial non-stationary approximating one. At the core of this al- 
gorithm there is a procedure which, for a given fixed non-interpolatory subdivision 
symbol a^'^^z), k > 0, effectively constructs a corresponding interpolatory symbol 
denoted by m^'^\z). The procedure is applied step-by-step for = 0, 1, . . .. For the 
sake of notational simplicity we can therefore omit the superscript k by denoting 
a^^^(z) = a{z) and m^^\z) = m{z). The construction stems from a theoretical result 
presented in [5l Theorem 2] which describes the conditions being satisfied for the 
associated interpolatory symbol m{z). In Subsection 13.11 this result is reviewed and 
generalized to some extent by removing unnecessary restrictions on the input symbol 
a{z). In the case where a{z) is of the form (12.51) and it is known in factorized form 
by means of the set of zeros {Oi, Te}e=i^... ^n, then an efficient method for computing 
a suitable representation of m{z) is described in Subsection 13. 2[ Finally, by putting 
all these ingredients together, Appint is formally stated in Subsection 13.31 

3.1 From approximating to interpolatory subdivision sym- 
bols 

In the matrix environment the linear operator Sa defined in (12. ip and associated with 
the symbol a{z) = '^i^^CLiz'', z & C \ {0} is represented by a bi-infinite Toeplitz- 
like matrix Sa = (ai_2j), i,j € Z. Since a{z) is a Laurent polynomial, say a{z) = 
'YTj=-K^j^'' 1 i^^lk-Kl) kftl} > 0) it follows that Sa is banded with bandwidth 
[|] at most. Let p{z) = Yl'j=-hPj^'^ ^ ^^^{\P-h\, \Ph\} > 0, be another Laurent 
polynomial and denote by V the bi-infinite Toeplitz matrix associated with p{z), 
namely, V = {pi~j)- Observe that V is again banded with bandwidth h. For the 
product operator 

S: = V ■ Sa = (sij), i,j e Z, 

we have 

i+h h 

= Pi-raj—2j = PlO,i-2'j-l = ■5j+2j+l, G Z. 

r=i—h £=—h 

This means that the product operator 5 is a bi-infinite Toeplitz-like matrix of the 
same form as the subdivision operator Sa with entries Sij = Si-2j, i,j G Z. By 
setting 

h+K 

q{z) = a{z) ■ p{z) = qjz^, {qj = ii \j\ > h + k), 

j=-h-K 

we find that 

h 

Qj =YPi o-j-i, -{h + k) < j < h + K, 

i=—h 
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and, therefore, 

Qi~2j = Si J = Sj_2j, i,j E Z. 

There follows that the product operator S can be seen as the subdivision operator 
associated with the Laurent polynomial q{z), i.e., 

S = q{z) = a{z) ■ p{z), 

where a{z) is the symbol of and p{z) can be suitably chosen in such a way to 
satisfy the interpolation condition. By expressing q{z) in terms of its sub-symbols 

q{z) = q even(^^) + ^ " ^ odd(^^) zeC\ {0}, 

we find that 

q{z) + q{-z) = 2-q even(^^)- 

Then by imposing the interpolation condition (12. 6p . i.e., q gyen^"^) = 1) arrive 
at the relation 

a{z) -piz) + a{-z) -pi-z) = 2 (3.8) 

which is a generalized Bezout equation providing necessary and sufficient conditions 
for a Laurent polynomial p{z) to convert the subdivision operator associated with 
a{z) into the interpolating subdivision operator generated by q{z) = a{z) -pi^z). 

Suitable coefficient-wise representations ofp{z) are introduced to investigate con- 
ditions under which the (generalized) Bezout equation is solvable as well as to de- 
velop effective computational methods for its solution. Observe that if p{z) is of the 
form 

p{z) = p^z^ + p.+,z^+' + ...+ p.+mZ^^"', (3.9) 
with m = 2k — 1, and, moreover, it satisfies 

a{z) ■ p{z) + {-iya{-z) ■ p{-z) = 2z\ < j < 2m, + 1, (3.10) 

then z~^p{^z^ solves (13.81) . Computing polynomial solutions of fl3.10p of the form 
(13. 9p reduces in a matrix setting to solving a structured hnear system whose co- 
efficient matrix is Sylvester-like. Let = [a_K, . . . , Oq, . . . , a^]^ G R^''"'"-^ denote 
the coefficient vector of the Laurent polynomial a{z). The associated extended co- 
efficient vector a+ G M^K+m+i jg (jggng^j i^j = [o-o, 0, . . . , O] . Similarly let us 
introduce the extended coefficient vector a_ G M^K+m+i aggociated with the polyno- 
mial a{—z). Moreover let Z = {zij) G ]R2(m+i)x2(m+i) -^^ ^j-^g down-shift matrix given 
by Zij = Si-i,j, where 6ij is the Kronecker delta symbol. Set TI+ G ]]j2(m+i)x{m+i) 
the striped Toeplitz matrix 

7^+ = [a+|Za+|...|Z™a+], 

and, similarly, define 

7^_ = [a_|Za_|...|Z™a_]. 

The coefficient matrix of the linear system fl3inl) is 7^+ = [7^+|7^_] G M2(m+i)x2{m+i) 
or TZ~ = [JZ+\ — TZ-] G K2(m+i)x2{m+i) depending on the parity of j. It is well known 
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that TZ'^ and TZ~ are resultant matrices and, therefore, they are invertible if and only 
if a{z) and a{—z) are relatively prime polynomials. 

Due to the special structure of the polynomial pair {a{z),a{—z)) it is shown 
that both linear systems can be reduced to smaller systems of half the size. Let 
Pm+i e K2(m+i)x2(m+i)^ P^+i = (<^j,cr(j)) be the permutatiou matrix associated with 
the "perfect shuffle" permutation given by 



a : {l,...,2m+2}^{l,...,2m+2}, a{j) 



(j + l)/2 + m + 1, if j is odd; 
j'/2, if j is even. 

Furthermore, let Gm+i ^ ]R2fcx2A: ^y^q matrix defined by 



G. 



m+l 



Im+1 


—Dm+1 


Drn+1 


I m+l 



where An+i = diag[-l, (-l)^, . . . , {-lf~^, There follows that 

Pm+l ■ TZ ■ Cm+l = ^ © "H, 

where "H G M(™-+^)^(™+^) is a certain matrix and 



(3.11) 





-K- 


hl 











-K- 


h3 




a 






-K- 


h5 


0--K+4 


a 


-K+3 







0-K+2m+l Ct-K+2m 0,-K+2m-l 



Similarly we find that 

P„+i ■ 7^+ ■ G^Vi = ^ © (3.12) 

where "H G is a certain matrix, © denotes the direct sum with respect 

to the main anti-diagonal, and, moreover. 



a-K 









fl-K + 2 













fl-K+2 





. 0'-K+2m 0'-K+2m-l 0'-K+2m-2 ■ ■ ■ 

In this way we arrive at the following generalization of O Theorem 2]. 

Proposition 1. Let a{z) = z'^a{z) be a degree-n polynomial, n = m + 1, relatively 
prime with a{—z). Then and are invertible and, moreover, the polynomial 
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Pi{z), -k G {+, — }, with coefficients given by the entries of the i-th column of (71*) ^, 
1 < i < n, is the unique polynomial of degree less than n such that 

a{z)pl{z) a{-z)p^i{-z) = 2z^'~^\ l<i<n, ^ G {+, -} (3.13) 

where 

' 2 if * = +; 
1, elsewhere. 

As an immediate consequence of Proposition [1] we obtain the following. 

Proposition 2. Given a degree-n polynomial a{z) relatively prime with a{—z) and 
such that a(l) = 2, a(— 1) = 0, then the Laurent polynomials 

?Wg£). !<,<„. (3.14) 

where p*{z) solves (I3.13P . * G {+, — }, are the associated interpolatory symbols and 
satisfy 

m*(l) = 2, m*(-l) = 0, l<i<n. 

Remark 1. It is worth noting that Proposition^ defines a double family of associated 
interpolatory symbols depending on the sign of -k. In the symmetric case where a{z) 
is a symmetric polynomial, that is, aj = a^j, < j < the number of associated 
symbols halves since all the matrices "H, "H, H"*" and "H" are suitably related and, in 
particular, T-C^ can be obtained from Ti' by reversion of rows and columns. 

These results provide a practical way to construct a family of finitely supported 
interpolatory masks from a given approximating one consisting in computing the 
matrix and reading its entries. This approach seems to be especially tai- 

lored for symmetric Hurwitz subdivision symbols which result into computations 
with totally positive (TP) Hurwitz matrices. The procedures described in [20] can 
be adjusted for the efficient and stable computations of the coefficients of the inter- 
polatory masks generated in the B-spline case and "shifted" affine combinations of 
them (see [5l Section 4]). However, in the case of exponential B-splines and their 
affine combinations the approximating symbol is generally known by assigning the 
spectrum of the symbol, that is, its zeros with their multiplicity. It is therefore 
interesting to design a completely different machinery for solving fl3.13p using the 
information on the roots. 

3.2 A root-based polynomial equation solver 

Let us suppose that 



m 



a[z) = oq + aiz + . . . +anz'^ = an]^(z - zj^ 

j=0 
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with Zi 7^ Zj ifi ^ j and ko + . . . + km = n. Then it is shown that the unique solution 
Pi{z) of (13.131) can be obtained by imposing certain interpolation conditions at the 
zeros of a(z) and a{—z). 

Let us start by recalling the concept of Hermite- Lagrange interpolation poly- 
nomial of a given differentiable function f{z) on the set of nodes rjo,. . . ,rii with 
multiplicities Hq, . . . ,hi, ho + . . . + hi = r + 1, respectively. Suppose that the func- 
tion f{z) possesses derivatives f^^\r]i), < j < hi — 1, < i < i. Then there 
exists a unique polynomial Hf{z) of degree at most r satisfying the interpolation 
conditions 

HfiVi) = f^'\m). < J < /i, - 1, < ^ < 

This polynomial is generally referred to as the Hermite-Lagrange interpolation poly- 
nomial of f{z) on the prescribed set of nodes. By setting uj{z) := {z — tiqY^ ■ ■ ■ [z — 
rji)^'^ we find the Lagrange-type representation 



and, equivalently, the partial-fraction representation 

i=0 s=l ^ \j=0 / i=0 s=l ^ 

where 

and, moreover, by Leibniz's rule 

i=0 •' ^ «V / / z=r\i 

Let £ = 2m -M and r/o = ^o, • • • , '7(«-i)/2 = ^m, ??(£+i)/2 = -^^o, •••,'?£ = -^^m with 
multiplicities /iq = ^(f+i)/2 = ko, . . . , hi = /i(£_i)/2 = A;^- Observe that 

r + 1 = /iQ . . . + /i£ = 2A;o + . . . + 2A;m = 2n 

and 

m m 

i=0 i=0 

By replacing the right hand side f{z) = 2z^^~^* of fl3.13p . where t is fixed and 
1 <t < n, with its Hermite-Lagrange form we find that 

"Vol-:) J tJ^rt(^-w)' 

Since 0(2;) and a{—z) are relatively prime we can separate the partial fraction de- 
compositions of the two rational functions on the left-hand side. This gives the 
following 
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Proposition 3. Let a{z) = anYYjLoi^ ~ ^jY^ ^ polynomial of degree n, where 
Zi ^ Zj if i ^ j , ko + . . . + km = n and a{z) and a{—z) are relatively prime. Then, 
the unique polynomial solution Pt{z), 1 < t < n, -k & {+, — }, of (13. 13^ satisfies 



j=0 

where 

1 / 2^2*-^* \^^^ 

oo{z) is the monic polynomial associated with a{z)a{—z) and Wi{z) = ^ ^ . 

Example 1. To illustrate the computational meaning of the previous result let us 
consider the interpolatory symbols associated with the cubic exponential B-spline with 
k-level symbol 

where the parameter v^^^ G (0, +oo) is defined through the expression 



2 

with 6 G {6e, ^ = 1, A^}, as in LemmaUl As shown in this means that B^^\z) 
corresponds to h2. 5|) with N = 3, 9i = 0, 02 = t, 6^ = —t and ti = 2, T2 = t-^ = 1, 
and, moreover, once assigned the starting value v^~^'^ G (— l,+oo), the parameter 
v^^'' can be recursively updated at each successive iteration through the formula 



^(fc)^ /!^!:!1±1^ A;>0. (3.15) 



For any fixed k > , the symmetric interpolatory scheme of smallest support associ- 
ated with B^\z) is obtained from the choice i = 2 and -k = — in (I3.13p . By using 
Proposition\^we find that the corresponding solution P2{z) is given by 



p(^\z) = \, , ' , — ^ = — (-z^ + 2 t;^ + l)z - 1) , 

which from PropositionlE defines the interpolatory symbol 

The partial fraction decomposition is not a flexible computational tool and sev- 
eral difliculties arise in order to flnd eflicient updating procedures for computing the 
solutions of (I3.13P associated with slightly modifled symbols (as usually it is the 
case in non-stationary subdivision schemes depending on a parameter, see Section 
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[5]). In this respect the tool of incomplete partial fraction decomposition [H] is much 
more suited. The general strategy proceeds as follows. From the partial fraction 
decomposition we get two polynomials h{z) and k{z) of degree less than n such that 

1 h(z) kiz) 
^ r = ^ + — r. Since aiz) is given in factored form we can determine 

a[z)a[—z) a[z) a[—z) 

k{z) as the Hermite-Lagrange polynomial interpolating the function g{z) = l/a{z) 
on the zeros of a{—z). Then the polynomial Pt{z) which solves f l3.13p can be ob- 
tained by means of the polynomial division between 2z^^~^ k{z) and a{—z). Again 
this operation reduces to computing the Hermite-Lagrange polynomial interpolat- 
ing 22;^*"^ ■ k{z) on the zeros oia{—z). In the case where the initial symbol a{z) is 
modified by a linear or a quadratic factor, both the two steps in the above procedure 
can be modified accordingly. For instance the polynomial k{z) can be specified in 
the form k{z) = ki{z) + a{—z)'ilj{z), where ki{z) is the Hermite-Lagrange polyno- 
mial interpolating the function g{z) = l/a{z) on the zeros of a{—z) and ip{z) is 
a linear factor whose coefficients are determined so that k{z) satisfies the modified 
equation. This approach has been implemented and used for computing the interpo- 
latory symbols associated with certain affine combinations of exponential B-splines. 
Some computational results are shown in Section [5l 



3.3 The Appint algorithm for the non- stationary case 

So far we have introduced a quite general strategy for deriving a family of interpola- 
tory symbols from a given approximating symbol based on the solution of equation 
f l3.13p . In the non-stationary setting, we compute a family of non- stationary interpo- 
latory subdivision schemes associated with a non-stationary approximating one via 
the solution of (13.131) at each recursion step. Therefore, the procedure we consider 
turns out to be as follows: assuming {a^''\z), k > 0} are the degTee-n{k) symbols of 
an approximating non-stationary scheme with a^'^^z) and a^''\—z) relatively prime 
for all k > OjWe construct the non- stationary interpolatory subdivision scheme based 
on the symbols {m[^l^{z), k > 0} where, for each k, m^^^l^{z), 1 < i{k) < n{k), is 
one of the interpolatory symbols satisfying (13.131) . Here and hereafter for the sake 
of simplicity we omit the superscript * G {+, — } since we assume that the sequence 
{i{k),'k), A; > 0, is given in input and, therefore, m['^l.^{z) denotes the unique solution 
of (I3.13P for the given pair {i{k),-k). Surely, the performance of the non-stationary 
subdivision scheme will depend on the selection of the sequence (z(/c),T<r), k > 0. The 
computational kernel consists of finding the solution of (13.131) for the input symbol 
a^''\z) and the fixed pair (z(A;),*). This task can be accomplished by the inversion 
of the corresponding matrices "H* or, alternatively, by means of the procedure de- 
scribed in the previous section based on computing the incomplete partial fraction 
decomposition. The auxiliary routine Solve takes in input a suitable representation 
of a^^^(2;) together with the pair (i{k),-k) and returns as output the corresponding so- 



lution pl^Liz) of (13.131) . For clarity we describe the overall procedure in algorithmic 



form. 
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Appint Algorithm 
Input: {a^''\z), k > 0}, degree-n(/c) symbols; 

{{i{k),ic), k > 0}, with 1 < i{k) < n{k) 
For A; = 0,1,... 

Check whether a^'^\z) is relatively prime with a^''\—z) 
Set pf^l^iz) :=Solve[aW(z), 

Construct the interpolatory symbol mi^l,^{z): = — ^^m- 
Output: {m'^hz), k>0} 



Some theoretical properties of the computed sequence {m\i^l-^{z), A; > 0} are 
discussed in Section H] whereas computational examples are reported in Section [51 



4 Properties of non-stationary interpolatory sub- 
division schemes derived from their approxi- 
mating counterparts 

For the family of non- stationary interpolatory subdivision schemes generated by 
symbols {ml {z),k > 0}, 1 < i < n{k), we can prove an important reproduction 
result: the exponential polynomial space reproduced by the interpolatory scheme 
is the same function space generated by the approximating scheme it is originated 
from. To prove it, we first need a preliminary result given in [12]. Within the rest 

of this section Vt,^ is the space given in Definition [1] and z\ := e~2fc+i , £ = 
1,---,A^, A;>0.' 



Proposition 4. Let {mS''\z), k>Q} he a sequence of interpolatory symbols. The 
subdivision scheme associated with such a sequence reproduces Vt,j if and only if 
for each k >0 

m(^)(zf^) = 2, m('^)(-^f^) = 0, i = l,---,N 

(4.16) 

|^m('=)(±^f ) =0, r = 1,--- ,r,-l, i=l,---,N. 
We are now in a position to state the reproduction result. 

Proposition 5. Let {a^''\z), k > 0} be a sequence of symbols with'a^^\z) relatively 
prime withaS'^\—z) for all k >0. If the non-stationary approximating subdivision 
scheme based on the symbols {a^^\z), k > 0} generates the space Vt^^, then for 
all 1 < i < n{k) the non-stationary interpolatory subdivision scheme based on the 
symbols 

whenever convergent, reproduces the same space Vr.^- 
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Proof: Due to [25l Theorem 1] the symbols a*-^^ (z) satisfy 

a(^)(-zf)=0, £^a('=)(-.f) = 0, r = l,---,r,-l, £=l,---,iV. 

By the Leibnitz's differentiation rule, we easily get an analogous relation to be 
satisfied by all mf'\z) (for any 1 < i < n{k)) that is 

mf)(-zf)=0, £-mf)(-zf ) = 0, r = l,---,r,-l, £=l,---,iV. 

It remains to consider the behavior of m[''\z) and its derivatives at the points zf^\ 
Now, since for each k 



it follows that 



as well as 



m\''\z) + mf =2, l<i< n{k), 
m\ '[z\ = 2 



= {-ir'^A\-^f^) = 0, r = 1, . . . , r, - 1, £ = 1, ■ ■ ■ , AT. 
The use of Proposition H] concludes the proof. □ 

Remark 2. We notice that, if an interpolatory subdivision scheme is Vp^^- gen crating, 
then due to the interpolatory nature ( that is due to the fulfillment of equation (2.1)), 
it is also Vt,^ -reproducing. 

Remark 3. Unfortunately, contrary to the result in Proposition a general re- 
sult concerning convergence and/ or smoothness of a non- stationary interpolatory 
subdivision scheme induced by a non- stationary approximating one is not available. 
However, in all specific examples discussed in Section\^ and many others we tested, 
convergence and smoothness analysis of the induced non-stationary interpolatory 
subdivision schemes is provided. From the examples we see that the smoothness 
order of the interpolatory scheme is the half of that of the approximating one it 
is originated from. This observation gives us a hint for a theoretical result to be 
investigated in future researches. 

5 Interpolatory exponential reproducing non-sta- 
tionary subdivision schemes 

Aim of this section is to show the application of our strategy to a family of approxi- 
mating schemes depending on free parameters. This leads to a parameter-dependent 
family of corresponding interpolatory schemes that can be used to design interesting 
new non- stationary interpolatory schemes. In particular, we show that by means 
of a five term affine combination of exponential B-splines, we can generate novel 
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smooth non-stationary interpolatory subdivision schemes possessing very interest- 
ing reproduction properties. 

Let us consider the interpolatory scheme based on the symbols ml^l (z) introduced 

in Example [H The approximating scheme with symbols {B^''\z), k > 0} was 
originally introduced in [19] where the authors also showed its capability generation 
of the function space V^^^ = {1, x, e*^, e~^^} (see also [21]). According to the results 
in Section HJ the associated interpolatory scheme turns out to be the 4-point 
interpolatory scheme reproducing the function space V^^^ = {1, x, e*^, e~^^} (see also 
[21]). The reproduction properties of this scheme can be improved by considering the 
family of approximating subdivision schemes given by a 5-term affine combination 
of B^^^ (z) of the form 

aW(2) = 5f (2) (a'^"^ + -I- (1 - 2aW - 2/3^)^2 + ^^z^ + a^'^h") 

= B''^^ iz) (a<'''^ + /3"°' + V(4a"°>+/3"'')"-4aW , ^(fc) 2(l-2g(^) -4af') ) ^2] 

3 ^ ■'y 2 / V l3(>'> + ^/(:laW+l3('')y^-4ai'') J 

where a^'^\f3^'^'^ E M are free parameters. By imposing the primality conditions 
for a^''\z),a^''\—z) it turns out that fl3.13p can be solved whenever a'-'^-* 7^ and 
^ {0, 1 - 2«(^), ^&^^lz±p±l}_ In the case a^"'^ = the equation can be 
degree-reduced in such a way that a polynomial solution can still be found. 

By applying the procedure described in Subsection 13.21 we have computed the 
polynomial p^'^^z) corresponding with the pair {i{k),-k) = (4, — ), k >0, and set 

By accurately choosing the free parameters a^''^ and (3^''\ we can obtain an inter- 
polatory scheme m^''\z) that improves the properties of the interpolatory scheme 
mi^l{z) associated with the combined symbol B^\z). Improvements can concern 
with its reproduction capabilities and/or its smoothness order. In particular: 

1. When a^^) = and (3^'^^ = i, a^'^^z) = ("+^\'^g+f;'^''"+^\ namely it is the 

exponential B-spline that generates Vq^-^ = {1, x, x^, x^, e*^, e~*^}. The symbol 
mS^\z) is the interpolatory 6-point scheme that reproduces the same space 
(as previously shown in [21]). 



2. When a^^'^ = and /J^'^) = then 

.(fe) {z + 1)2(^2 + 2t;W^ + l){z^ + 2(2(t;('=))2 - l)z + 1) 



a^'>iz) 



16(wW)2(wW + 1) 



namely it is the exponential B-spline that generates Vq,^ = {1, x, e*^, e~*^, 
e^*^, e~2*^'}, while m^'^^z) is the interpolatory 6-point scheme that repro- 
duces the same space (see, again, [2T]). 



When a^^^ = and 13'-''^ = ^(j^^, then 



(,) ^ {z+inz' + 2v('^)z + iy 
^' 8(w(^) + l)2 
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tx 



namely it is the exponential B-spline that generates Vq^^ = e*^,e 
xe*^, xe~*^} and m^''\z) is the interpolatory 6-point scheme that reproduces 
the same space [21]. 



4. When = 8(,(.))2(,w+i)(2.(.)_,)2 and = 4(Sl?)?(2'w!i)' > then 

^(fe), . _ (z+l)^(z2+2^;(fe)z+l)(z^+2(4(t;(fe))3-3^(fe))z+l)(z^+2(2(^(fe))2-l)z+l) 
« [^) - 32(i;('=))2(v(fe) + l)2(2t;(fe)-l)2 



i.e. it is the exponential B-spline generating Vs,^ = {1, x, e*^, e ^'^j e^*^, e" 



-2ix 



^3tx^ e ^*^}, while m^'^\z) defines the interpolatory 8-point scheme that re- 
produces the same space (see [B]). 

5. When a^^^ = ^(^^(k)/^^^^^) and /3^^'^ = |^^^, we deal with the exponential 
B-spline 

^(,) , _iz + 1)2(^2 + 2v('h + 1)2(^2 + 2(2(t;W)2 - l)z + 1) 



32(t;W)2(^;W + 1)^ 



generating the function space V8,7 = {1, x, e*^, e~*^, e^*^, e"^*^^ xe*^, xe~*^}. 
The symbols m^''\z) define a interpolatory 8-point scheme that repro- 
duces the same space (see Proposition [5]) . The smoothness of the subdivision 
scheme {m^'''^{z), k > 0} can be obtained through asymptotical equivalence 
[9] with the Dubuc-Deslauriers 8-point interpolatory scheme [H [13] . 



The last non-stationary interpolatory subdivision scheme corresponds to a new 
proposal never presented in the literature. Other interesting proposals can be ob- 
tained by assigning different suitable values to the free parameters a'-'^^ and l3^''\ In 
all these kinds of interpolatory schemes^ by making the parameter v^~^^ local, namely 
by assuming a different parameter in correspondence of each edge qiQi+i of the 
starting polyline, we can combine the two important issues of local shape control 
and special functions reproduction. This means that, in the same limit curve, we can 
include an alternation of exponential polynomial pieces in those regions where the 
starting samples belong to one of these curves and smooth limit segments with local 
tension otherwise. Also, due to the recurrence relation (13.151) . the shape parameter 
v^^'^ turns out to be independent of the parametric values t^^\ thus reducing compu- 
tational costs of the algorithm. In addition to the general reasons discussed in the 
introduction, these properties contribute to make these interpolatory subdivision 
schemes more convenient with respect to the corresponding classical interpolatory 
methods. 



6 Conclusions and future work 

A novel approach has been presented for the computation of a family of interpo- 
latory non-stationary subdivision schemes from a non-stationary, non-interpolatory 
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one. The approach reduces the updating problem either to the inversion of certain 
structured matrices (which can be of Hurwitz type or Sylvester resultant matrices) 
or to the solution of certain Bezout-like polynomial equations. If the approximat- 
ing symbols are defined in terms of spectral information it is shown that the partial 
fraction decomposition provides an effective tool for solving these equations by yield- 
ing a representation of the associated interpolatory symbols in terms of generalized 
interpolating conditions. The newly constructed interpolatory schemes are capa- 
ble of reproducing the same exponential polynomial space as the one generated by 
the original approximating scheme. Although a general result concerning the rela- 
tionship between convergence and/or smoothness orders of the approximating and 
interpolatory schemes is not yet available, ad hoc techniques can be used by show- 
ing that in many cases the proposed approach leads to novel smooth non- stationary 
interpolatory subdivision schemes possessing very interesting reproduction proper- 
ties. The analysis of more general convergence properties of the subdivision schemes 
generated by our techniques is an ongoing research. 
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